/*
**  File: AquaAssign.c
**
**  This file is part of the AQUA program suite, which is being
**  developed in the NMR Spectroscopy Department, Bijvoet Center for 
**  Biomolecular Research, Utrecht, the Netherlands and the BioMagResBank, 
**  University of Wisconsin-Madison, USA.

**
**  VERSION   : 3.2
**  DATE      : January 2001
 
    MERGED INTO AQUA 1.5 on Dec 11, 1997

    AUTHOR(S) : Jurgen F. Doreleijers based on programs by Ton Rullmann
    PURPOSE   : see files doc/assign_method.txt and doc/assign_output.txt

    Notes:

    Redundant restraints can be extracted
    
    INSERTION CODES are neglected, but
    CHAIN identifiers are handled correctly.
    
    The restraints should not contain 
    * dubbels! (Checked by other aqua procedures)
    
    Restraints between atoms in one prochiral (stereoassignable) group will
    be omitted from analysis.
    Aug 23, 1996 and UPDATE May 29, 1998 (also when writing the restraints).
    
    Restraints including pseudo-atoms for protons which are also listed
    seperately in the list no longer cause a problem.
    Nov 5, 1996

*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include "Compiler.h"
#include "General.h"

#include "AquaMacros.h"
#include "AquaTypes.h"
#include "AquaTypes_Assign.h"

/****************  function prototypes  ****************/

#include "AquaFuncts.h"
#include "AquaFuncts_dr.h"
#include "AquaStrucset.h"
#include "AquaFiles.h"
#include "MenuFuncts.h"
#include "Functs.h"

#include "AquaAssign.h"
#include "AquaCommon.h"
#include "nrutil.h"


/****************  globals  ****************/

extern Boolean  iaf, batch;                           /* import from Functs */
extern int      terseness;                            /* import from Functs */
char            *program_name = "AquaAssign";         /* export to Functs */
char            *program_vers = aqua_version;         /* export to Functs */
static Card     line, descr, inp;
char            temp_value[4];

/* My own globals */
static Mol		*molptr;				    	    
static Atom		*atoms;
static int		struc_count;
static Distinfl		*dist_infl;
static Distrst		*dist_rst;
static Distvall		*dist_vall; 
static int		dist_count, num_atoms;
static Extra_Distinfl	*extra_distinfl;           /* Types defined in: AquaTypes_Assign.h */
static Triplet		triplet[MAX_TRIPLETS], ext_triplet[MAX_TRIPLETS];
static int		tri_count, ext_tri_count;

enum {	opt_ext, opt_ass, opt_readass, opt_sort, opt_energy };

/* BEGIN_EXTH */

int main( int argc, char *argv[] )
{
    char         *fct_name = "AquaAssign main";
    Strucset     theSet;
    Model        *models = NULL;
    DistdataList inf_list, val_list;
    DistinfTyp   ftype = inf_extended;
    Qtyp         qtype, qtv;
    int          s_count, count;
    Filnam       dinf_filnam, drst_filnam, dval_filnam;
    Filnam       filnam_s, filnam_i, filnam_r;
    Card         ident;
    int          item;

/* My own variables */
    int	      todo_count, ass_count;
    
    enum { item_do };
    static MenuItem menulist[] =
    {
        { item_EXEC,  "exec",  none, on },
        { item_SET,   "set",   none, on },
        { item_do,    "do",    none, on }, 
        { item_QUIT,  "quit",  none, on } 
    };
    static OptionItem optionlist[] =
    {
        { opt_batch,       "batch",       logical, on }, 
        { opt_terseness,   "terseness",   integer, on }, 
        { opt_prtaq,       "print",       logical, on }, 
        { opt_ext,         "useext",      logical, on }, 
        { opt_ass,	   "ass",         logical, on }, 
        { opt_readass,     "readass",     logical, on }, 
        { opt_sort,        "sort",        logical, on },
        { opt_multi,       "multiple",    logical, on },
        { opt_prtlib,      "printlib",    logical, on }, 
        { opt_prtres,      "printres",    logical, on }, 
        { opt_guess,       "guess",       logical, on }, 
        { opt_energy,      "energy",      logical, on }
    };

    static Menu menu =
    {
        "Main menu", "command", "Give command:", 
        sizeof(menulist)/sizeof(MenuItem)
    };
    static Menu options =
    {
        "Options menu", "option", "Give option:", 
        sizeof(optionlist)/sizeof(MenuItem)
    }; 


/* initialize */

    ProgramInit( argc, argv );
    ProjectInit();
    StrucsetInit( &theSet, "theSet", true );
    SetMenu( &menu, menulist );
    SetOptions( &options, optionlist );
    ModOption( opt_ass  , "on"  );
    ModOption( opt_multi, "on" );
    sprintf( temp_value, "%d", terseness );
    ModOption( opt_terseness, temp_value );    

    inf_list.count = 0;
    val_list.count = 0;

/* the main menu loop */

    while ( true )
    {
        item = ReadMenuItem();
        switch ( item )
        {
            case item_EXEC:
		DoExec();
                break;
            case item_SET:
		ReadOption();
		batch = OptionTrue( opt_batch );
                if ( !Check_General_Options () ) {
	            printf( "Error in main: general options set are not okay\n" );
	            MError( );
                }
                terseness   = GetOptionIVal ( opt_terseness );                
                break;
            case item_do:
		goto doit;
            case item_QUIT:
            case item_NULL:
                exit( 0 );
                break;
            default:
                ProcessMenuError( menu_fatal );
                break;
        }
    }

/* do the main thing */

doit:

/* process the distance-values file saving the header lines and file names */

    *dval_filnam = '\0';
    ProcessDistvalFile( &val_list, &dist_count, &struc_count, &models, &qtv, 
                        theSet.filnam, dinf_filnam, drst_filnam,
		        ident, true, &ftype, dval_filnam );
    if ( ftype != inf_extended )
    {
	FATAL( " * Wrong file/data type - Code error - case 1" );
    }
    dist_vall = (Distvall *) (val_list.lst_desc->data);

/* process the distance-restraints file (file name read from distance-values file) */

    if ( IsEmptyStr( drst_filnam ) )
    {
	nrerror( " * Name of dist-restr file not found on distance-values file");
    }
    dist_rst = ProcessDistrstFile( &count, &qtype, filnam_s, filnam_i, 
                                   ident, false, drst_filnam );
    if ( NOTEQ( theSet.filnam, filnam_s ) ||
         NOTEQ( dinf_filnam, filnam_i )   ||
         dist_count != count )
    {
        FATAL(" * Inconsistent information on values and restraints files");
    }

/* process the distance-info file (file name read from distance-values file) */

    ProcessDistinfFile( &inf_list, &count, &s_count, &qtype, filnam_s, filnam_r,
                        ident, false, &ftype, dinf_filnam );
    if ( ftype != inf_extended )
    {
	FATAL( " * Wrong file/data type - Code error - case 2" );
    }
    dist_infl = (Distinfl *) (inf_list.lst_desc->related->data);
    if ( NOTEQ( theSet.filnam, filnam_s ) ||
         NOTEQ( dinf_filnam, filnam_i )   ||
         dist_count != count )
    {
        nrerror( " * Inconsistent information on values and info files" );
    }

/* open and read the strucset file */
    if ( terseness < 6 ) {
        printf("Opened and read the strucset file\n");
    }

    if ( EMPTYSET( theSet ) )
    {
	if ( !ProcessStrucsetFile( &theSet, -1 ) )
	    exit( 1 );
    }
    
    molptr	= theSet.molecule;
    atoms	= molptr->atoms;
    num_atoms	= molptr->atom_count;
    
    
    
/* Assign specific routines */

    if ( dist_count > MAX_RESTRAINTS )
	nrerror("Error too many restraints :-)");
    if ( struc_count > MAX_MODELS )
	nrerror("Error too many models :-)");
    
    Get_Triplets();

    if ( OptionTrue( opt_readass ) )
	Read_AssignDAT();

    Check_Triplets( &todo_count );

    Init_Extra_Distinfl();

    Assign_Triplets( todo_count, &ass_count );
    Sort_Triplets( tri_count );
    if ( terseness < 9 ) {
        Print_Triplets( todo_count, ass_count );  
    }
       
    Write_Corrected_Restraints();

    Write_Assignments(ass_count);

    if ( terseness < 6 ) {
        printf("Finished AquaAssign\n\n");
    }
    exit(0);
}


static void Get_Triplets(void)
{
/* 
 * Get the assignable pairs (triplets) & involved restraints
 * This will check if	- an atom is in a restraint
 *			- the atom is assignable according to dist_infl
 * The counter partner will always be initialized too.
 */
    Distidx  *C;
    int i;
    int tri_idx, count_both, countC, idxC;
    FILE    *fp = stdout;
      
    if ( terseness < 6 ) {
        printf("Looking for triplets\n");
    }
    
    tri_count = 0;
    
    for ( i=0; i < dist_count; i++ ) {	 

	if ( dist_infl[i].A[0].idx == dist_infl[i].B[0].idx ) {
	    
            if ( terseness < 9 ) {
	        printf("Restraint between: ");
	        Write_Atom( fp, atoms + dist_infl[i].A[0].idx );
	        printf(" and its stereo-pair or itself will be neglected.\n");
            }
/*
	    nrerror("Error in Get_triplets: see restrictions on restraints");                                     
*/
	    continue; /* Do not make a triplet for it */
	}
	
	count_both = 2;
	while ( count_both-- > 0) {
	
	    if (count_both == 1) {
		countC = dist_infl[i].count_A;
		idxC   = dist_infl[i].A[0].idx; /*Only first element considered */
		C      = dist_infl[i].A;  /* HB2 should always be before HB3  */
	    } else {
		countC = dist_infl[i].count_B;
		idxC   = dist_infl[i].B[0].idx;
		C      = dist_infl[i].B;		
	    }
	    
	    
	    if ( ( countC < 1 ) || ( countC > MAX_ATOMS_TRIPLET ) )
		nrerror( "ERROR wrong number of atoms in list" );
	    
	    if ( countC > 1 ) {         /* Its perhaps possible to assign it  */
		
		if ( !Atom_In_A_Triplet(idxC, &tri_idx) ) {      
							       /* New triplet */
		    if (++tri_count > MAX_TRIPLETS) 
			nrerror( "ERROR too many assignable triplets" );

		    triplet[tri_count-1].P		    = C[0].idx;
		    triplet[tri_count-1].Q		    = C[1].idx;
		    triplet[tri_count-1].X		    = C[2].idx;
		    triplet[tri_count-1].assigned	    = false;
		    triplet[tri_count-1].ext_assigned	    = false;
		    triplet[tri_count-1].swapped	    = false;
		    triplet[tri_count-1].restrt_index[0]    = i;
		    triplet[tri_count-1].restrt_count	    = 1;
		    
		} else {                     /* Add as an old triplet member  */

		    if ( ++(triplet[tri_idx].restrt_count) > 
			MAX_RESTRAINTS_PER_TRIPLET ) 
			nrerror( "ERROR too many restraints for this triplet" );
			
		    triplet[tri_idx].restrt_index    [ triplet[tri_idx].restrt_count -1 ] = i;
		}
	    }
	}
    }
    if ( terseness < 6 ) {
        printf("Found %d triplets\n", tri_count);
    }
} 

static Boolean Atom_In_A_Triplet (int idx_atom, int *tri_idx) 
{                       /* Returns the triplet index in which the atom occurs */

    Boolean found = false;
        
    *tri_idx = -1;
    
    while ( ( !found ) && ( ++(*tri_idx) < tri_count ) ) 
	if ( triplet[*tri_idx].P == idx_atom )
	    found = true;
	    	
    return( found );
}


static void Check_Triplets ( int *todo_count )
{
    int i, j;
    Boolean swapped;
    
/* Checks if triplets are already assigned
 * 
 * throws out restraints which are between members of one triplet eg Gly HA2/HA3
 * 
 * Sorts accordingly to (1) todo or not and (2) to residue/name order.
 */
    
    if ( terseness < 6 ) {
        printf("Start checking %d triplets\n", tri_count);
    }
    *todo_count = 0;
    
    if ( tri_count > MAX_TRIPLETS )
	nrerror( "ERROR in Check_Triplets (-1)" );
	
    for (i=0;i<tri_count;i++) {
    
	if ( (triplet[i].P > num_atoms) || (triplet[i].Q > num_atoms) ||
	     (triplet[i].X > num_atoms) )
	    nrerror( "ERROR in Check_Triplets (0)" );
	    
	if ( (triplet[i].restrt_count < 1 ) ||
	     (triplet[i].restrt_count > MAX_RESTRAINTS_PER_TRIPLET ) )
	    nrerror( "ERROR in Check_Triplets (1)" );
	    
	for (j=0;j<triplet[i].restrt_count;j++) {
	    	    
	    if ( ( triplet[i].restrt_index[j] < 0 ) || 
	         ( triplet[i].restrt_index[j] >= dist_count  ) )
		nrerror( "ERROR in Check_Triplets (2)" );


	    if ( ( dist_infl[ (triplet[i].restrt_index[j]) ].count_A == 3 ) &&
	         ( dist_infl[ (triplet[i].restrt_index[j]) ].count_B == 3 ) )
		if ( dist_infl[ (triplet[i].restrt_index[j]) ].A[2].idx ==
		     dist_infl[ (triplet[i].restrt_index[j]) ].B[2].idx ) 
		    nrerror( "ERROR in Check_Triplets (3)" );
	}
	
	if ( External_Assigned( triplet[i].P, &swapped ) ) {
	    triplet[i].ext_assigned = true;
	    triplet[i].assigned	    = true;
	    triplet[i].swapped	    = swapped;
	    triplet[i].ass_idx	    = -1;    /* Flag showing when assigned */
	} else {
	    triplet[i].ass_idx	    = 0;    
	    (*todo_count)++;
	}
	
    }
    Partition_Triplets (*todo_count);
    Sort_Triplets      (*todo_count);

    if ( terseness < 6 ) {
        printf("Will try to assign %d triplets\n", *todo_count);
        printf("Finish checking routine\n");
    }
    
}

static Boolean External_Assigned(int atom_idx, Boolean *swapped)
{     

/*  Check against info from assignment file 

    Returns true if the atom with the given atom index is in a
    externally assigned triplet. If so it also returns the swapped 
    status 
*/
    
    int i;
    *swapped = false;
    
    if ( ( atom_idx < 0 ) || ( atom_idx >= num_atoms ) )
	nrerror( "ERROR in routine External_Assigned; bad atom_idx" );
	
    for (i=0;i<ext_tri_count;i++) {
	if ( atom_idx == ext_triplet[i].P ) 
	    if ( ext_triplet[i].ext_assigned == true ) {
		*swapped = ext_triplet[i].swapped;
		return(true);
	    } else
		return(false);
    }
    
    return(false);
}

static void Partition_Triplets(int todo_count)
{                /* orders the not external ass. triplets first in the row  */
    int i, j, last_ass, count;
    Boolean found;
    
    if ( terseness < 6 ) {
        printf("Start partitioning the triplets accordingly\n");
        printf("%d (out of %d) triplets to really do; putting them first\n", 
	        todo_count, tri_count);
        if (todo_count==0) 
	     printf("Will do nothing\n");
    }

    i	     = 0;
    count    = 0;
    last_ass = tri_count;
    
    while ( (i < last_ass ) && (count < todo_count) ) {
	if (triplet[i].ext_assigned == true) {    /* Needs to be moved to end */
	    j = last_ass;
	    found = false;
	    while (!found) {
		j--;
		if (j==i) 
		    nrerror( "ERROR (0) in routine Partition_Triplets" );
		if ( triplet[j].ext_assigned == false ) {
                    if ( terseness < 6 ) {
		        printf("Swapping triplets %i and %i\n", i, j);
                    }
		    Swap_Triplets(i, j);
		    found = true;
		    last_ass = j;
		    count++;	    
		}
	    }
	} else {
	    count++;
	}
	if (++i > tri_count ) 
	    nrerror( "ERROR (1) in routine Partition_Triplets" );
    }
    if ( terseness < 6 ) {
        printf("Finish partition routine\n");
    }
}

static void Swap_Triplets(int i, int j)
{
    Triplet tt;
    
    if ( (i==j) || (i<0) || (j<0) || (i >= tri_count) || (j >= tri_count) )
	nrerror( "ERROR in routine Swap_Triplets" );
	
    tt = triplet[i];
    triplet[i] = triplet[j];
    triplet[j] = tt;
}

static void Sort_Triplets(int todo_count)
{       

/* Sorts first "todo_count" number of triplets according to 
    chain ID, residue number and atom index
    using the "piksrt" algorithm of Numerical Recipes.
    Yes, I know it's slow (N^2).
*/
	    
    Triplet temp;
    int	    i, j;
    
    if ( terseness < 6 ) {
        printf("Start sorting %d triplets\n", todo_count);
    }
    if ( tri_count < todo_count ){
	nrerror( "ERROR(s) in sorting the triplets, too many ;-)" );
    }
    
    for (j=1; j<todo_count; j++) {
	temp = triplet[j];
	i=j-1;
	while ( (i>=0) && Triplet_Order(temp, triplet[i]) ) {
	    triplet[i+1] = triplet[i];
	    i--;
	}
	triplet[i+1] = temp;
    }
    if ( terseness < 6 ) {
        printf("Finish sorting\n");
    }
}

static Boolean Triplet_Order (Triplet tri_1, Triplet tri_2)
{
     /* Returns true if first atom in triplet 1 should come before that of 2 */

    int	resnum1, resnum2, comp;
    
    if ( tri_1.P == tri_2.P )
	nrerror( "ERROR in Triplet_Order" );

    resnum1 = (atoms + tri_1.P)->res_num_orig;
    resnum2 = (atoms + tri_2.P)->res_num_orig;

    comp = memcmp ( (atoms + tri_1.P)->chain, 
		    (atoms + tri_2.P)->chain, size_ID );
    
    if ( comp < 0 )
	return(true);
    if ( (comp == 0) && (resnum1 <  resnum2) )
	return(true);
    if ( (comp == 0) && (resnum1 == resnum2) && (tri_1.P < tri_2.P) )
	return(true);
		
    return(false);
}

static void Print_Triplets (int todo_count, int ass_count)
{
/* Output look-alike:
 In general it is a table containing a fixed number of colums. The last given
 assignment state is the same as the internal one.

PS	pseudo Atom representing the triplet
RES	residue
#	number of ..
Ex	external assignment state
In	internal '		    '
Initial first round
Final	second round
pM	% Models favoring the state
pE	average (over all models and restraints) relative 
	(with respect to unfavorable state, negated) energy difference
dE	average energy difference (in Angstrom^2)
A       assignment state
I	insertion code
C	chain identifier

*/

    FILE    *fp = stdout;
    char    kar[4];
    int	    i, tot_swapped;
    double  hi_energy, lo_energy, dif_energy;
    double  per_energy, per_ass, per_swap, per_model, sum_hi=0.0, sum_lo=0.0;
    Atom    *at;
    
    printf("Assignment statistics:\n");

    printf("Number of triplets              :  %2d\n",  tri_count);
    printf("Number to assign                :  %2d\n", todo_count);
    printf("Number assigned                 :  %2d\n",  ass_count);
    
    if ( todo_count != 0 )
	per_ass = (100.0 *  ass_count / (double) todo_count );
    else
    	per_ass = 0;
	
    
    printf("Percentage (assigned/to assign) : %3.0f\n", per_ass);
    
    
    tot_swapped = 0;
    for (i=0;i<tri_count;i++)
	if ( triplet[i].ext_assigned || triplet[i].assigned )
	    if ( triplet[i].swapped )
		tot_swapped++;

    printf("Number swapped                  :  %2d\n", tot_swapped );
    
    if ( todo_count != 0 )
	per_swap = (100.0 *  tot_swapped / (double) tri_count );
    else
    	per_swap = 0;
	
    printf("Percentage swapped              : %3.0f\n", per_swap);

    if ( terseness < 6 ) {
        printf("\nAssignment table:\n");

        printf("Triplet:                            *                                     *\n");
        printf("                         Assigned   * Statistics                          *\n");
        printf("            Chain        | Swapped  *           Energies                  *\n");
        printf("RES     #   C    PS      # Ex+In |  * pM  pE    delta   high    low  A  # *\n");
        printf("***************************************************************************\n");

        for (i=0;i<tri_count;i++) {

	    at = atoms + triplet[i].X;

	    Write_Atom( fp, at ); 

	    kar[0] = Boolean_To_Char(triplet[i].ext_assigned);
	    kar[1] = Boolean_To_Char(triplet[i].assigned);

	    printf(" %3d %c  %c  ", triplet[i].ass_idx, kar[0], kar[1]);

	    if ( triplet[i].ext_assigned || triplet[i].assigned )
	        kar[0] = Boolean_To_Char(triplet[i].swapped);
	    else
	        kar[0] = '-';

	    printf("%c   ", kar[0] );

	    if (triplet[i].ext_assigned)
		    printf("\n");
	    else {

	        per_model =  (100 * (double) triplet[i].int_num_mdl_fav ) / 
				         (double) struc_count;

	        if ( (per_model < 0) || (per_model > 100) ) {
		    printf ( "ERROR in Print_Triplets (wrong per_model = %d)\n",  (int) per_model);
		    nrerror( "ERROR in Print_Triplets (wrong per_model)" );		
	        }

	        if ( triplet[i].int_energy[0] > triplet[i].int_energy[1] ) {
		    hi_energy = triplet[i].int_energy[0];
		    lo_energy = triplet[i].int_energy[1];
	        }
	        else {
		    hi_energy = triplet[i].int_energy[1];
		    lo_energy = triplet[i].int_energy[0];
	        }

	        if ( hi_energy < SMALL_POS_DOUBLE_NUMBER ) 
		     hi_energy = SMALL_POS_DOUBLE_NUMBER;

	        dif_energy = hi_energy - lo_energy;

	        per_energy = 100 * (dif_energy / hi_energy);

	        if ( ( per_energy < 0 ) || ( per_energy > 100 ) )
		    nrerror( "ERROR in Print_Triplets (wrong per_energy)" );


	        kar[0] = Boolean_To_Char(triplet[i].int_ass);

	        printf("%3.0f %3.0f %7.3f %7.3f %7.3f %c %2d\n", per_model, per_energy, 
		        dif_energy, hi_energy, lo_energy, kar[0], triplet[i].restrt_count);

	        sum_hi += hi_energy;		    
	        sum_lo += lo_energy;		    
	    }
        }

        printf("***************************************************************************\n");
        printf("Violation high states : %8.3f (A^2/model) \n", sum_hi/struc_count);
        printf("Violation low states  : %8.3f (A^2/model) \n", sum_lo/struc_count);
    }
}

  
static char Boolean_To_Char(Boolean boo)
{
    
    if (boo) 
	return('T'); 
    else 
	if (boo != false) {
	    nrerror( "ERROR in Boolean_To_Char" );
	    return('N'); /* Will never be executed */
	} else
	    return('F');
}

static Boolean Char_To_Boolean(char ch)
{
    ch = toupper(ch);
    if (ch == 'T') 
	return(true);
    if (ch == 'F') 
	return(false);
	
    nrerror( "ERROR in Char_To_Boolean" );
    
    return(false); /* Will never be executed */
}


static void Write_Assignments (int ass_count)
{

    FILE    *fp;
    int	    i;
    char    PathName[] = ".";
    char    File[]="AssignDAT-new";
    Filnam  FileName;
    char    ch;
    Atom    *at1, *at2;
    Qtyp    qtype = assign_data;
    
    sprintf(FileName,"%s/%s", PathName, File);
    F_OPEN_PROTECTED (fp,FileName,"w");

    if ( terseness < 6 ) {
        printf("Writing Assignment file %s\n", FileName);
    }

    WriteAquaIDLines( fp, qtype, FileName );

    fprintf(fp, "$ REMARK: Format used for the AQUA Assignment Table.\n");
    fprintf(fp, "$ REMARK: \n");
    fprintf(fp, "$ REMARK: Authors: JF Doreleijers and JAC Rullmann\n");
    fprintf(fp, "$ REMARK: AQUA names based on:\n");
    fprintf(fp, "$ REMARK: \"Recommendations for presentation of NMR structures of proteins\n");
    fprintf(fp, "$ REMARK:  and nucleic acids\"\n");
    fprintf(fp, "$ REMARK: IUPAC-IUBMB-IUPAB Interunion Task Group (draft version 11/29/95)\n");
    fprintf(fp, "$ REMARK: \n");
        
    fprintf(fp, "$ COUNT: %d\n", tri_count);
    fprintf(fp, "$ ASSIGNED_COUNT: %d\n", ass_count + ext_tri_count);
    fprintf(fp, "$ BEGIN_DATA:\n");

    for (i=0;i<tri_count;i++) {
    
	at1     =  atoms + triplet[i].P ;
	at2     =  atoms + triplet[i].Q ;
	
	if (triplet[i].assigned) 
	    ch = Boolean_To_Char(!triplet[i].swapped);
	else 
	    ch = '?';
	
	fprintf (fp, "def  ");
	Write_Atom(fp, at1);
	fprintf (fp, "  ");
	FPRINTID(fp, at2->atom_int ); 
	fprintf (fp, "     %c\n", ch);

    }
    
    fclose(fp);
}


static void Assign_Triplets(int todo_count, int *ass_count)
{               /* Assigns the triplets where possible. See AquaTypes_Assign.h */

    Boolean continue_ass, *tri_done;
    int	    pick, num_loops, max_loops, old_loop_ass_count;
    
    max_loops = 1;
    num_loops = 0;
    *ass_count = 0;

    tri_done = CREATE_NEW(Boolean, todo_count);       /* Parallel with triplet array */
    
    if ( OptionTrue( opt_ass ) )
	continue_ass = true;
    else
	continue_ass = false;
    
    while ( continue_ass && ((*ass_count) < todo_count) ) {
    
	if (++num_loops > max_loops)
	    nrerror( "ERROR in Assign_Triplets; too many loops " );
	    
        if ( terseness < 6 ) {
	    printf("\n **** Doing ASSIGN LOOP: %d *************\n\n", num_loops);
        }
	continue_ass = false;
	old_loop_ass_count = *ass_count;
	Init_Boolean_Array ( tri_done, todo_count, false ); 
	                            /* Used for processing status of TRIPLET  */

	while ( Pick_Triplet( tri_done, &pick, todo_count ) ) {

	    if ( Assign_A_Triplet( pick, ass_count ) ) {
		continue_ass = true;
                if ( terseness < 6 ) {
		    printf("ASSIGNED ! (%3d = total)\n", *ass_count);
                }
	    }
	    
	    tri_done[ pick ] = true;

	}
		
        if ( terseness < 6 ) {
	    printf("\nFOUND %3d ASSIGNMENTS IN LOOP %d.\n", 
		*ass_count - old_loop_ass_count, num_loops);
        }

    }
    
    if ( terseness < 9 ) {
        printf("Number of assignments: %2d\n", *ass_count);
    }
    if ( terseness < 6 ) {
        printf("Finished assignment routine\n\n");
    }
}


static Boolean Pick_Triplet(Boolean *tri_done, int *best, int todo_count)
{

/* Looks for the 'best' triplet in the list 
 * This is a triplet which has not been done (picked) or assigned before.
 * Returns true and the index of the best triplet and 
 * returns false and -1 when no triplets remain.
 */
 
    int i, temp, contacts, contacts_assigned;
    FILE *fp=stdout;
    
    contacts = 0;
    contacts_assigned = 0;
    *best = -1;

    for ( i=0; i<todo_count; i++ ) {
    
	if ( !( tri_done[i] || triplet[i].assigned ) ) {
    	    
	    temp = Number_Assigned_Contacts( i );
	    
	    if ( triplet[i].restrt_count > contacts ) {
	    
		*best = i;	
		contacts = triplet[i].restrt_count;
		contacts_assigned = temp;
	    } else    
		if ( ( triplet[i].restrt_count == contacts ) &&
		     ( temp > contacts_assigned ) )            {
		     
		    *best = i;	
		    contacts = triplet[i].restrt_count;
		    contacts_assigned = temp;
		}
	}
    }
    
    if ( *best == -1 )
	return ( false );

	
    if ( terseness < 6 ) {
        printf( "Picked triplet %3d ( %3d contacts, %3d assigned ) ", 
	    *best, contacts, contacts_assigned );
        Write_Atom(fp, atoms + triplet[*best].X); 
        printf("\n");
    }
	    
    return ( true );
}

static int Number_Assigned_Contacts(int tri_idx)
{

/*  Will return the number of contacts of the triplet in which 
    the other partner (not in this triplet) is assigned 
    First a check will be made if the first or second atom occurs in triplet
*/

    int i, rstrnt_idx, assigned_count;
    Boolean A_has_triplet, B_has_triplet;
    
    assigned_count = 0;
    
    if ( ( tri_idx >= tri_count ) || ( tri_idx < 0 ) )
	nrerror( "ERROR in routine Number_Assigned_Contacts (0)" );
        
    for ( i=0; i <   triplet[tri_idx].restrt_count; i++ ) {
    	
	rstrnt_idx = triplet[tri_idx].restrt_index[i];
	
	if ( dist_infl[rstrnt_idx].A[0].idx == triplet[tri_idx].P )
	    A_has_triplet = true;
	else
	    A_has_triplet = false;

	if ( dist_infl[rstrnt_idx].B[0].idx == triplet[tri_idx].P )
	    B_has_triplet = true;
	else
	    B_has_triplet = false;
	    
	/* Do you remember this:
	 * 
	 * !(A&&B) = (!A)||(!B) en
	 * !(A||B) = (!A)&&(!B)
	 * 
	 * Can be checked by the waarheidstabllen like:
	 * 
	 *     A
	 *     t f
	 * B  tF T  if the operation is !(A&&B)
	 *    fT T 
	 * 
	 */
	 
	if ( A_has_triplet && B_has_triplet ) {
	    printf("Useless restraint found in triplet %d\n", tri_idx);
	    nrerror( "ERROR in routine Number_Assigned_Contacts (2)" );
	}
	if ( ! (A_has_triplet || B_has_triplet) )
	    nrerror( "ERROR in routine Number_Assigned_Contacts (3)" );


	if ( A_has_triplet && extra_distinfl[rstrnt_idx].boo_B )
	    assigned_count++;
	    
	if ( B_has_triplet && extra_distinfl[rstrnt_idx].boo_A )
	    assigned_count++;

    }
       
    return ( assigned_count );
}

static void Init_Extra_Distinfl(void)
{

/*  Executed after reading ext. assignments but before doing assignments.
    Fills each extra_distinfl with its boo_'s and ass's and cor's 
    An atom in an unassigned triplet should have assignment F and 2 */

    int	    i, tri_idx_A, tri_idx_B;
    Boolean Assigned_A, Assigned_B;
    
    extra_distinfl = CREATE_NEW( Extra_Distinfl, dist_count );      /* Only now used */    

    for (i=0;i<dist_count;i++) {
    
	Assigned_A = Atom_Assigned( dist_infl[i].A[0].idx, &tri_idx_A );
	Assigned_B = Atom_Assigned( dist_infl[i].B[0].idx, &tri_idx_B );
		
		
	if ( Assigned_A )				    /* A */
	    extra_distinfl[i].boo_A = true;
	else
	    extra_distinfl[i].boo_A = false;
	    	    
	if ( tri_idx_A == -1 )
	    extra_distinfl[i].ass_A = 0;
	else 
	
	    if ( Assigned_A )
	    
		if ( triplet[tri_idx_A].swapped )  
		    extra_distinfl[i].ass_A = Swapped_Value(dist_infl[i].ass_A);
		else     
		    extra_distinfl[i].ass_A = dist_infl[i].ass_A;
		    
	    else
		extra_distinfl[i].ass_A = 2;
	
	
	if ( Assigned_B )				    /* B */
	    extra_distinfl[i].boo_B = true;
	else
	    extra_distinfl[i].boo_B = false;
	    	    
	if ( tri_idx_B == -1 )
	    extra_distinfl[i].ass_B = 0;
	else 
	
	    if ( Assigned_B )
	    
		if ( triplet[tri_idx_B].swapped )  
		    extra_distinfl[i].ass_B = Swapped_Value(dist_infl[i].ass_B);
		else     
		    extra_distinfl[i].ass_B = dist_infl[i].ass_B;
	    else
		extra_distinfl[i].ass_B = 2;
		

    }
    Calculate_Corrections ( extra_distinfl, dist_infl, dist_count);    
}

                
static Boolean Atom_Assigned(int at_idx, int *tri_idx)
{            

/* This routine will return two things:
    (1) If the atom is assigned it returns true otherwise false.
    (2) If it is contained in a triplet its index will be returned otherwise -1.
*/

    int i;
    
    for (i=0;i<tri_count;i++) {
    
	if (triplet[i].Q == at_idx)                   /* Checks second atom */
	    nrerror( "ERROR in routine Atom_Assigned" ); 
	
	if (triplet[i].P == at_idx) {              /* Only checks first atom */
	    *tri_idx = i;
	    if ( triplet[i].assigned )
		return (true);
	    else
		return (false);
	}    
	
    }
    
    *tri_idx = -1;
    return (true);                                 /* not found in a triplet */
}


static void Init_Boolean_Array (Boolean *array, int count, Boolean value)
{
    int i;
    
    for ( i=0; i<count; i++ ) 
	array[i] = value;
}

static Boolean Assign_A_Triplet (int pick, int *ass_count) 
       /* Checks if the given triplet can be assigned and does so if possible */   
{
    int i, j, rstrt_idx;
    double high_energy, diff_energy;
    Boolean swap_state, energy_criterium, model_criterium;
    double  **per_mdl_st;

    per_mdl_st = dmatrix( 0, struc_count-1, 0, 1 );
    Fill_Dmatrix (per_mdl_st, 0, struc_count-1, 0, 1, 0.0);
    
    for (i=0;i<MAX_STATES;i++) {                          /* Loop over states */

	if (i==0)
	    swap_state = false;
	else
	    swap_state = true;
	
	triplet[pick].int_energy[i] = 0;
	
	for (j=0;j<triplet[pick].restrt_count;j++) {  /* Loop over restraints */
	    rstrt_idx = triplet[pick].restrt_index[j];
	    triplet[pick].int_energy[i] +=
		Energy_Restraint( rstrt_idx, pick, swap_state, per_mdl_st);
	}
	
    }
    
    triplet[pick].int_num_mdl_fav = 0;
    
    for (i=0;i<struc_count;i++) {              /* Decision on ini_num_mdl_fav */
	if (per_mdl_st[i][0] <= per_mdl_st[i][1])
	    triplet[pick].int_num_mdl_fav++;
		
	if ( OptionTrue( opt_energy ) )
	    printf("Energy (%3d %4d) = %8.3f, %8.3f\n", 
		pick, i, per_mdl_st[i][0], per_mdl_st[i][1]);		 
    }
    
    diff_energy = triplet[pick].int_energy[0] - triplet[pick].int_energy[1];

    if (diff_energy > 0) {                       /* Decision on favored state */
	triplet[pick].int_fav_swap = true;                         /* Swapped */
	triplet[pick].int_num_mdl_fav = struc_count 
					- triplet[pick].int_num_mdl_fav;
	high_energy = triplet[pick].int_energy[0];
    } else {
	triplet[pick].int_fav_swap = false;
	high_energy = triplet[pick].int_energy[1];
    }
    
    if ( high_energy == 0) high_energy = 0.001;
    if ( struc_count == 0)
	nrerror("ERROR in routine Assign_A_Triplet: No structures");

    
    if ( ( ( my_fabs(diff_energy) / high_energy) >= ENERGY_REL_CRITERIUM ) &&
           ( my_fabs(diff_energy) >= ENERGY_ABS_CRITERIUM * struc_count ) )
	energy_criterium = true;
    else
	energy_criterium = false;

    if ( ( triplet[pick].int_num_mdl_fav / 
         ( (double) struc_count ) ) >= MODEL_CRITERIUM )
	model_criterium = true ;
    else
    	model_criterium = false ;

    if ( true ) { /* always assign to unswapped state unless criteria fullfilled */
    
	(*ass_count)++;                                        
	
	/* Adjust the triplet */
	
	triplet[pick].assigned = true;
	triplet[pick].int_ass  = true;
	triplet[pick].ass_idx  = *ass_count;

	if ( energy_criterium && model_criterium && 
           ( triplet[pick].int_fav_swap == true) ) 
	    triplet[pick].swapped  = true;
	else
	    triplet[pick].swapped  = false;


	for (j=0;j<triplet[pick].restrt_count;j++) {  /* Loop over restraints */
	
	    rstrt_idx = triplet[pick].restrt_index[j];

	    if ( triplet[pick].P == dist_infl[rstrt_idx].A[0].idx )
	    {				     
		extra_distinfl[rstrt_idx].boo_A = true;
		
		if ( triplet[pick].swapped )
		    extra_distinfl[rstrt_idx].ass_A = 
				       Swapped_Value (dist_infl[rstrt_idx].ass_A);
		else 
		    extra_distinfl[rstrt_idx].ass_A = dist_infl[rstrt_idx].ass_A;
				    
	    } else {
		
		if ( triplet[pick].P != dist_infl[rstrt_idx].B[0].idx )
		    nrerror( "ERROR in Assign_A_Triplet; wrong triplet info" );
		    
		extra_distinfl[rstrt_idx].boo_B = true;
		
		if ( triplet[pick].swapped )
		    extra_distinfl[rstrt_idx].ass_B = 
				        Swapped_Value (dist_infl[rstrt_idx].ass_B);
		else 
		     extra_distinfl[rstrt_idx].ass_B = dist_infl[rstrt_idx].ass_B;
	    }
	}
    }
    
    free_dmatrix( per_mdl_st, 0, struc_count - 1, 0, 1 );

    return (triplet[pick].assigned);
}

static double Energy_Restraint (int rstrt_idx, int triplet_idx, 
		                Boolean swap_state, double **per_mdl_st)
{
	/* Calculates the energy contribution of one restraint over all models
		    given the assignment state of the triplet to be assigned. */
		    
    double  energy, energy_mdl, upper, lower;
    int	    i, values_idx, idxA, idxB, dummy, state;
    double  correction = 0;

    if (rstrt_idx >= dist_count) 
	nrerror( "ERROR in routine Energy_Restraint; too large restrt_index" );
    
    /* Find values index belonging to the assignment states of both */
    /* First a few checks */
    
    if ( extra_distinfl[rstrt_idx].boo_A && 
	 extra_distinfl[rstrt_idx].boo_B    )
	     nrerror( "ERROR in Energy_Restraint; restraint totally assigned" );
	     
    
    if ( dist_vall[rstrt_idx].count != (dist_infl[rstrt_idx].count_A * 
					dist_infl[rstrt_idx].count_B ) )
	nrerror( "ERROR in routine Energy_restrt; unmatched counts " );
			  
    if ( !Triplet_Restraint(rstrt_idx, triplet_idx, &dummy) )
	nrerror("ERROR in Energy_Restraint; triplet doesn't contain restrt");
	
    if ( triplet[triplet_idx].P == dist_infl[rstrt_idx].A[0].idx ) {

	if ( triplet[triplet_idx].P == dist_infl[rstrt_idx].B[0].idx )
	    nrerror("ERROR in Energy_Restraint; useless restraint found");

	if ( swap_state )
	    idxA = Swapped_Value( dist_infl[rstrt_idx].ass_A );
	else
	    idxA = dist_infl[rstrt_idx].ass_A;
	
	idxB = extra_distinfl[ rstrt_idx ].ass_B;
	
	if (idxB == 2) {
	    correction = extra_distinfl[rstrt_idx].cor_B;
	    if (correction > (BAD_CORRECT - SMALL_POS_DOUBLE_NUMBER) ) {
		nrerror("ERROR in Energy_Restraint; bad corr B");
	    }
	}

    } else {
    
	if ( swap_state ) 
	    idxB = Swapped_Value( dist_infl[rstrt_idx].ass_B );
	else
	    idxB = dist_infl[rstrt_idx].ass_B;
	
	idxA = extra_distinfl[ rstrt_idx ].ass_A;
	
	if (idxA == 2) {
	    correction = extra_distinfl[rstrt_idx].cor_A;
	    if (correction > BAD_CORRECT - SMALL_POS_DOUBLE_NUMBER )
		nrerror("ERROR in Energy_Restraint; bad corr A");
	}
    }
    
    values_idx = idxA * dist_infl[rstrt_idx].count_B + idxB;
    

    upper =  dist_rst[rstrt_idx].up_bnd + correction;    
    lower =  dist_rst[rstrt_idx].lo_bnd - correction;
    if (lower < SMALL_POS_DOUBLE_NUMBER ) 
	lower = SMALL_POS_DOUBLE_NUMBER ;
    
    energy = 0;
    
    if (swap_state)
	state = 1;
    else
	state = 0;
    
    for (i=0;i<struc_count;i++) {
    
	energy_mdl = Energy_Model(i, rstrt_idx, values_idx, upper, lower);

	energy += energy_mdl;
	per_mdl_st[i][state] += energy_mdl;
	
    }
    
    return (energy);
}

static double Energy_Model ( int struct_idx, int rstrt_idx, int val_idx, 
			     double upper, double lower )
{
    double value, diff;
    int idx;
    
    idx = struct_idx * dist_count + rstrt_idx;
    value =  dist_vall[idx].values[val_idx];
    
    if ( dist_rst[rstrt_idx].typ_bnd == 'b' || dist_rst[rstrt_idx].typ_bnd == 'B' )
	if (value > upper) {
	    diff = value - upper;
	    return( diff * diff );
	} else
	    if (value < lower) {
		diff = lower - value;
		return( diff * diff );
	    } else 
		return( 0 );   
    
    if ( dist_rst[rstrt_idx].typ_bnd == 'u' ) 
	if (value > upper) {
	    diff = value - upper;
	    return( diff * diff );
	} else 
	    return( 0 );
    
    if ( dist_rst[rstrt_idx].typ_bnd == 'l' ) 
	if (value < lower) {
	    diff = lower - value;
	    return( diff * diff );
	} else 
	    return( 0 );
    
    nrerror( "ERROR in routine Energy_Model; wrong type of bound" );
    return(-999.9);					    /* Never executed */
}
	
	
static int Swapped_Value (int value)
{
    switch ( value ) {
	case 0 : return (1);
	case 1 : return (0);
	default : nrerror( "ERROR in routine Swapped_Value" );
    }
    
    return (9); /* Will never be executed */
}


static Boolean Triplet_Restraint(int rstrt_idx, int tri_idx, int *idx)
{      
/* Returns true and the restraint index (as located in triplet) if in triplet 
   returns false and a restraint index of -1 if not present in triplet */
	    
    int i;
    	    
    *idx = -1;
    
    for ( i=0; i<triplet[tri_idx].restrt_count; i++)
	if ( triplet[tri_idx].restrt_index[i] == rstrt_idx ) {
	    *idx = i;
	    return ( true );
	}
	
    return ( false );  /* Restraint doesn't contain a triplet member */ 
}


static void Fill_Dmatrix (double **per_mdl_st, int i1, int i2, int j1, int j2, double x)
{
    int i, j;
    
    for (i=i1;i<=i2;i++)
	for (j=j1;j<=j2;j++)
	    per_mdl_st[i][j] = x;
}

static void Write_Atom( FILE *fp, Atom *at)
{
    /* Write atom in specific format */
    
    FPRINTID(fp, at->res_nam);
    
    fprintf(fp, " %3d   ", at->res_num_orig);
    
    FPRINTID(fp, at->chain);
    
    FPRINTID(fp, at->atom_int);
}



static void Read_AssignDAT(void)
{

/****************  Read_AssignDAT  ****************/
/*
**  Process the assignment definitions from Aqua's AssignDAT file.
**
**  Format assumed:
**  (1) standard Aqua header
**  (2) standard Aqua header
**  (3) $ COUNT: <triplet_count>
**  (4) $ ASSIGNED_COUNT: <triplet_count>
**  (5) $ BEGIN_DATA:
**  (6) def <res_name> <res_num> <insert-code> <chain_code> 
**	<atom_name_1> <atom_name_2> <assign_flag>
**  (repeat 6)
**
**  Not assigned triplets will be discarded.
**  Lines starting with # are treated as comments.
**  Empty or blank lines force end of reading.
**  Other lines generate an error.
**
**  NOTE:   res_num is treated as res_num_orig
**	    insertion codes can not be read.
Two types of input are supported; without 
def  GLN     2        HB2    HB3       T
or with chain ident:
def  ASN     8   A    HD21   HD22      ?
Wildcards are not allowed.
**
*/
    FILE     *fp;
    int      i, index;
    int	     triplet_count,  assigned_count, not_assigned_count;
    int	     triplet_count_read,  assigned_count_read;
    Boolean  done = false;
    Card     txt[6];
    Qtyp     qtype;
    int	    resnum;
    char    flag;
    Filnam  filnam;
    ID	    atnam1, resnam, chain;
    
    qtype = assign_data;
    ext_tri_count = 0;
    
    triplet_count = 0;
    assigned_count = 0;
    not_assigned_count = 0;
    
    fp = OpenAquaFile( qtype, "r", filnam, false );
    if ( !fp )
	nrerror("ERROR: in Read_AssignDAT not opened");
    if ( !ReadAquaIDLines( fp ) )
	nrerror("ERROR: in Read_AssignDAT no AquaIDLines");
        
    printf("Reading assignments\n");
    
/* process file header */

    while ( !done )
    {
	if ( ! ProcessAquaLine( line, descr, inp, size_LINE, fp ) )
	{
	    if ( !iaf )
		nrerror("ERROR: in Read_AssignDAT unknown");
	    nrerror("ERROR: in Read_AssignDAT empty file");
	}
        if      ( EQn( descr, "COUNT", 5 ) )
            sscanf( inp, "%i", &triplet_count_read );
        if      ( EQn( descr, "ASSIGNED_COUNT", 14 ) )
            sscanf( inp, "%i", &assigned_count_read );
        else if ( EQn( descr, "BEGIN_DATA", 10 ) )
            done = true;
    }

/* read the list */

    done = false;
    while ( !done )
    {
	if ( !ReadWord( line, size_LINE, true, fp ) )     /* input word; eol test */
	    break;
	if ( EQn( line, "#", 1 ) )
	{
	    done = FlushLine( fp ) < 0;
	    continue;
	}
        else if ( EQn( "def", line, 3 ) )                 /* start of def */
        {
	                                                  /* store atom def */
	    Blank( resnam,    size_ID );
	    Blank( atnam1,    size_ID );
	    Blank(  chain,    size_ID );

	    i = 0;
	  
	    while ( ReadWord( txt[i], size_LINE, true, fp ) )
		i++;
	    
	    switch (i) {
	    
		    case 5 : /* No chain ident */
			    /* printf(" A case 5 line \n"); */
			    CopyChars( resnam, txt[0], size_ID );
			    sscanf   (         txt[1], "%i", &resnum);
			    CopyChars( atnam1, txt[2], size_ID );
			    sscanf   (         txt[4], "%c", &flag);
			    break;
		    case 6 :
			    /* printf(" A case 6 line \n"); */
			    CopyChars( resnam, txt[0], size_ID );
			    sscanf   (         txt[1], "%i", &resnum);
			    CopyChars( chain,  txt[2], size_ID );
			    CopyChars( atnam1, txt[3], size_ID );
			    sscanf   (         txt[5], "%c", &flag);
			    break;
		    default:
			if (i < 5)
			    nrerror("ERROR: in Read_AssignDAT not enough words");
			if (i > 6)
			    nrerror("ERROR: in Read_AssignDAT too many words");
	    }
	    
	    triplet_count++;
	    
	    flag = toupper(flag);

	    if ( flag == '?' )
		    not_assigned_count++;
	    else {
		if (flag=='T' || flag=='F') {
		
		    assigned_count++;
		    
		    if ( ++ext_tri_count > MAX_TRIPLETS )
			nrerror("ERROR: in Read_AssignDAT (0) too many triplets" );
			
		    ext_triplet[ext_tri_count-1].ext_assigned = true;
		    ext_triplet[ext_tri_count-1].swapped = !Char_To_Boolean(flag);
		    
		    index = Atom_Index(resnum, atnam1, chain);
		    
		    ext_triplet[ext_tri_count-1].P = index;
		    ext_triplet[ext_tri_count-1].Q = -1; /* These aren't used */
		    ext_triplet[ext_tri_count-1].X = -1;
		    
		} else
		    nrerror("ERROR: in Read_AssignDAT, bad assignment flag"); 
	    }
	    
	    done = FlushLine( fp ) < 0;
        }
	else
	{
	    puts( "Unexpected line:" );
	    puts( line );
	    nrerror("ERROR: in Read_AssignDAT ");
	}
    }
    
    if ( triplet_count != triplet_count_read )
	nrerror("ERROR: in Read_AssignDAT, wrong counts (0)");

    if ( assigned_count != assigned_count_read )
	nrerror("ERROR: in Read_AssignDAT, wrong counts (0)");
	
    if ( not_assigned_count != triplet_count - assigned_count)
	nrerror("ERROR: in Read_AssignDAT, wrong counts (1)");

    fp = stdout;

    if ( terseness < 6 ) {
        printf("Finished reading assignment file\n");
        printf("Read %i assignments\n", assigned_count);
        
        for (i=0;i<assigned_count;i++) {

	    printf("Triplet %3d with first atom: ", i);

	    Write_Atom( fp, atoms + ext_triplet[i].P );

	    printf(" with flag: %c\n", Boolean_To_Char( !ext_triplet[i].swapped ) );
        }
        printf("Finished Read_AssignDAT procedure\n");
    }
        
}


static int Atom_Index(int resnum, ID atonam, ID chain)
{

/* Finds the atom index given its residue number, atom name and chain identifier.
   All should be internal                                                     */

    int	    i;
    Atom    *at;
    	
    at = atoms;
        
    for ( i=0; i<num_atoms; i++ ) {
	    
	if ( ( at->res_num_orig == resnum ) && 
	     ( EQID( atonam, at->atom_int) ) &&
	     ( EQID( chain,  at->chain) ) )
    	    return (i);
	    
	at++;
    }
    
    printf("Residue number: %3d\n", resnum);
    printf("Atom name "); PRINTID(atonam);
    printf(" and Chain "); PRINTID(chain); printf("\n");
    
    nrerror("ERROR: in AtomIdx,  couldn't find it");
    return(-1);                                          /* Never executed */
}


static void Write_Corrected_Restraints (void)
{

/* Write the corrected restraints according to the assignment and options set */

#ifndef MAX_SAMES
#define MAX_SAMES    4
#endif

    int	    i, j, rstrnt_count, *list, list_count;
    int     exit_flag, multi_restr_count;
    
    DRList          DR_list;
    DRestraint	    *restraint;
    Distinfl	    *infptr;
    Extra_Distinfl  *extra_infptr;
    Distrst	    *rstptr;
    Boolean	    *done, boo_A, boo_B;
    double	    correction;
    
    if ( terseness < 6 ) {
        printf("Start writing corrected restraints\n");
    }
    
    restraint =  CREATE_NEW(DRestraint, dist_count);
    InitDRestraints ( restraint, dist_count);
    
    DR_list.list = restraint;
    DR_list.struc_count = 1;

    done = CREATE_NEW(Boolean, dist_count);
    list = CREATE_NEW(int, MAX_SAMES);
    
    rstrnt_count = 0;
    extra_infptr = extra_distinfl;
    infptr       = dist_infl;
    rstptr       = dist_rst;
    
    
    Init_Boolean_Array ( done, dist_count, false );
        
    for (i=0;i<dist_count;i++) {
    	
	if ( done[i] ) {
	    extra_infptr++;
	    infptr++;
	    rstptr++;	    
	    continue;
	}
		
	switch ( rstptr->typ_bnd ) {
	    case 'B' :    restraint->stat = bnd_uld; break;
	    case 'b' :    restraint->stat = bnd_uplo; break;
	    case 'l' :    restraint->stat = bnd_lower; break;
	    case 'u' :    restraint->stat = bnd_upper; break;
	    default : nrerror("Wrong restraint status of typ_bnd");
	}
			
	correction = 0;
	
	boo_A = extra_infptr->boo_A;
	boo_B = extra_infptr->boo_B;

	if (boo_A) {                                                     /* A */
	    restraint->idx_A = infptr->A[ extra_infptr->ass_A ].idx;
	    restraint->atom_A = *(atoms + restraint->idx_A );
	} else {
	    restraint->idx_A = infptr->A[2].idx;
	    restraint->atom_A = *(atoms + restraint->idx_A );
	    correction += extra_infptr->cor_A;
	}

	if (boo_B) {                                                     /* B */
	    restraint->idx_B = infptr->B[ extra_infptr->ass_B ].idx;
	    restraint->atom_B = *(atoms + restraint->idx_B );
	} else {
	    restraint->idx_B = infptr->B[2].idx;
	    restraint->atom_B = *(atoms + restraint->idx_B );
	    correction += extra_infptr->cor_B;
	}

	/* Do not store restraints that are between the same atoms 
	    UPDATE at May 29, 1998
	*/
 	
	if ( restraint->idx_A == restraint->idx_B ) {
	    extra_infptr++;
	    infptr++;
	    rstptr++;	    
	    continue;
	}
		

	restraint->lo_bnd = rstptr->lo_bnd;
	restraint->up_bnd = rstptr->up_bnd;
	restraint->dist = rstptr->dist;
/*
 Deleting these statements because they give problems to bounds being
 not recognized later on by e.g. AquaFuncts_xplor.c writting out the
 constraints        
	if ( restraint->lo_bnd <= TWO_PROTON_RADII + SMALL_POS_DOUBLE_NUMBER ) 
		restraint->lo_bnd = NO_LOWER_BOUND_TOKEN;
*/
        
	/* Look for "SAMES" which are restraints further on in the list
	 * which are to be collapsed together with this one due to lack of
	 * assignments. The total of sames can be from 1 to MAX_SAMES. 
	 * A list will be used to keep track. */
	 	
	
	list_count = 1;
	list[0] = i;
	
	if ( !(boo_A && boo_B) )
	    Find_Sames (list, done, i, boo_A, boo_B, &list_count);
	
	for (j=1;j<list_count;j++) {     /* Already initialized the first one */
	    done[ list[j] ] = true;
	    restraint->up_bnd = Best_Value(restraint->up_bnd, dist_rst[ list[j] ].up_bnd, false);
	    restraint->lo_bnd = Best_Value(restraint->lo_bnd, dist_rst[ list[j] ].lo_bnd, true);
	}
		
	restraint->up_bnd +=  correction;     /* Correction can be 0 */
	restraint->lo_bnd -=  correction;     /* depending on option */
	
/* see above
#
#	if ( restraint->lo_bnd <= TWO_PROTON_RADII + SMALL_POS_DOUBLE_NUMBER ) 
#	    restraint->lo_bnd = NO_LOWER_BOUND_TOKEN;
*/			
	restraint++;
	rstrnt_count++;
	extra_infptr++;
	infptr++;
	rstptr++;	
    }
    
    DR_list.dist_count = rstrnt_count;
    if ( terseness < 6 ) {
        printf("We are left with %i from the initial %i restraints\n", 
	    rstrnt_count, dist_count);
    }
    
    if ( OptionTrue( opt_sort ) )
	SortRestraints( &DR_list, true );

    multi_restr_count = CleanupDRList( &DR_list, 
            0, 
            &exit_flag );
  
    
    
    if ( terseness < 6 ) {
        printf("We are left with %i from the initial %i restraints\n", 
	        rstrnt_count - multi_restr_count, dist_count);
    }
                
    ProcessDROutput( &DR_list, 
            multi_restr_count, 
            "noe", false, molptr, struc_count);    

}

static double Best_Value(double val1, double val2, Boolean swap)
{
    double highest, lowest;
    
    if ( val1 > val2 ) {
	highest = val1;
	lowest  = val2;
    } else {
	highest = val2;
	lowest  = val1;
    }
    
    if ( swap )                          /* For lowerbound the lowest is best */
	Swap_Doubles(&highest, &lowest);

    return(lowest);
}

static void Swap_Doubles ( double *a, double *b )
{
    double x;
    
    x = *a;
    *a = *b;
    *b = x;
}

static double my_fabs ( double a )
{
    if ( a < 0 ) {
	return ( -a );
    } else {
	return (  a );
    }
}


static void Find_Sames(int *list, Boolean *done, int rst_idx, 
		Boolean boo_A, Boolean boo_B, int *list_count)
{

#define MAX_BOOLEANS 4

#ifndef MAX_SAMES
#define MAX_SAMES    4
#endif

    /* Finds the 'Same restraints' those only differing in triplet atoms */


    int	    i, j;
    Boolean found;
    FILE    *fp=stdout;

    
    for (i=rst_idx+1;i<dist_count;i++) {
    
	if ( done[i] )
	    continue;

	found = false;
	
	if ( !boo_A && !boo_B )
	
	    if ( ( dist_infl[i].count_A == 3 ) && ( dist_infl[i].count_B == 3 ) )
		if (( ( dist_infl[rst_idx].A[0].idx == 
		        dist_infl[     i ].A[0].idx ) &&
		      ( dist_infl[rst_idx].B[0].idx == 
		        dist_infl[     i ].B[0].idx ) ) || 
		    ( ( dist_infl[rst_idx].A[0].idx == 
		        dist_infl[     i ].B[0].idx ) &&
		      ( dist_infl[rst_idx].B[0].idx == 
		        dist_infl[     i ].A[0].idx ) ) )
		    found = true;  


	if ( !boo_A && boo_B ) {
	
	    if ( dist_infl[i].count_A == 3 )
		if ( ( dist_infl[rst_idx].A[                        0 ].idx == 
		       dist_infl[     i ].A[                        0 ].idx ) &&
		     ( dist_infl[rst_idx].B[ dist_infl[rst_idx].ass_B ].idx == 
		       dist_infl[     i ].B[ dist_infl[      i].ass_B ].idx ) )
		    found = true;

	    if ( dist_infl[i].count_B == 3 )
		if ( ( dist_infl[rst_idx].A[                        0 ].idx == 
		       dist_infl[     i ].B[                        0 ].idx ) &&
		     ( dist_infl[rst_idx].B[ dist_infl[rst_idx].ass_B ].idx == 
		       dist_infl[      i].A[ dist_infl[      i].ass_A ].idx ) )
		    found = true;
	}
	    
	    
	if ( boo_A && !boo_B ) {
	
	    if ( dist_infl[i].count_B == 3 )
		if ( ( dist_infl[rst_idx].B[                        0 ].idx == 
		       dist_infl[     i ].B[                        0 ].idx ) &&
		     ( dist_infl[rst_idx].A[ dist_infl[rst_idx].ass_A ].idx == 
		       dist_infl[     i ].A[ dist_infl[      i].ass_A ].idx ) )
		    found = true;
	    
	    if ( dist_infl[i].count_A == 3 )
		if ( ( dist_infl[rst_idx].B[                        0 ].idx == 
		       dist_infl[     i ].A[                        0 ].idx ) &&
		     ( dist_infl[rst_idx].A[ dist_infl[rst_idx].ass_A ].idx == 
		       dist_infl[     i ].B[ dist_infl[      i].ass_B ].idx ) )
		    found = true;

	    
	}

	if ( found ) {                                           /* Its the same .. */

	    if ( ++(*list_count) > MAX_SAMES )
		nrerror("ERROR in FindSame: too large list of sames");
		
	    list[*list_count - 1] = i;
	    
	}
    }
    	
    if ( ( (*list_count >         2) && ( boo_A ||  boo_B) ) ||
         ( (*list_count > MAX_SAMES) && (!boo_A && !boo_B) ) ) {

	printf("Wrong list:  \n");
	for (j=0;j<(*list_count);j++) {
	
	    printf("%4d ", list[j]);
	    Write_Atom ( fp, atoms + dist_infl[ list[j] ].A[0].idx );
	    printf("  ");
	    Write_Atom ( fp, atoms + dist_infl[ list[j] ].B[0].idx );
	    printf("\n");
	    
	}

	nrerror("ERROR in FindSame: got wrong list");
    }	 
}


/* END_EXTH ----------------------- only private functions follow */

/* Not used for now:
static void Show_Extra_Distinfl(void)
{

    Show the extra distance information
    int	    i;
    
    printf("Showing the extra distance information\n");
    
    printf("#  # Ass # # Correction\n");
    
    for (i=0;i<dist_count;i++)
	printf("%5d %c %c %1d %1d %8.3f %8.3f\n", i, 
	    Boolean_To_Char(extra_distinfl[i].boo_A), 
	    Boolean_To_Char(extra_distinfl[i].boo_B),
	    extra_distinfl[i].ass_A, extra_distinfl[i].ass_B, 
	    extra_distinfl[i].cor_A, extra_distinfl[i].cor_B   );    
}

static void Show_Distinfl(void)
{

    Show the distance information
    int	    i, idxA, idxB;
    FILE *fp=stdout;

    printf("Showing the distance information\n");
        
    for (i=0;i<dist_count;i++) {
    
	idxA = dist_infl[i].ass_A;
	idxB = dist_infl[i].ass_B;

	printf("%5d %1d %1d ", i, idxA, idxB );  
	
	Write_Atom( fp,  atoms + dist_infl[i].A[idxA].idx );
	
	printf(" ");
	
	Write_Atom( fp,  atoms + dist_infl[i].B[idxB].idx );
	
	printf("\n");
    }  
}

static void Show_Triplets(void)
{

    int i, j;
    FILE *fp=stdout;
    
    printf("Will show %d triplets \n", tri_count);
    printf("  # Atom X                Ass/Swap         ");
    printf("P    Q    X    #      restrt_count members\n"); 

    for ( i=0; i < tri_count; i++ ) {
    
	printf("%3d ", i);
	
	Write_Atom( fp,  atoms + triplet[i].X);
	
	printf(" %c %c %c %3d  ", 
	    Boolean_To_Char(triplet[i].assigned), 
	    Boolean_To_Char(triplet[i].ext_assigned), 
	    Boolean_To_Char(triplet[i].swapped), 
	    triplet[i].ass_idx);
	    
	printf(" %5d%5d%5d   %2d  ", 
	    triplet[i].P, triplet[i].Q, triplet[i].X, triplet[i].restrt_count);
	    
	for ( j=0; j < triplet[i].restrt_count; j++ ) {
	
	    if ( ( j%10 == 0) && ( j > 0 ) )
		printf("\n%61s", "");
		
	    printf("%5d", triplet[i].restrt_index[j]);
	}
	
	printf("\n");
    }    
}
*/
